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摘要 : 快速 传 里 叶 变 换 比 侍 里 叶 变 换 有 更 好 的 算法 性 能 ， 是 射电 干涉 成 像 基础 算法 , 但 
因为 天 线 阵列 的 不 规则 采样 , 需 使 用 网 格 化 算法 将 可 见 度数 据 重 采样 到 规则 的 网 格 上 , 才能 
应 用 快速 传 里 叶 变 换 。 基 于 卷 积 的 网 格 化 计算 具有 计算 密集 型 和 和 迭代 型 的 特点 ， 特 别 是 在 处 
理 海量 可 见 度数 据 的 情况 下 , 高 性 能 的 网 格 化 计算 对 加 速 整个 成 像 过 程 就 显得 尤为 重要 。 为 
了 缓解 数据 处 理 的 压力 , 在 现 有 处 理 整 块 数据 和 支持 多 核 计算 的 算法 基础 上 , 拓展 应 用 Dask 
并 行 计算 框架 ,不 仅 将 数据 分 块 并 分 配 到 多 线程 上 ， 提 高 数值 计算 效率 , 而且 动 态 的 分 布 式 
任务 调度 策略 优化 了 网 格 化 的 实时 处 理 。 实 验 结果 表明 多 核 CPU 利用 率 显 著 提高 ,即使 增加 
数据 量 , 也 能 进一步 提升 网 格 化 算法 性 能 。 分布 式 任务 调度 能 够 将 ( 单 ) 多 测量 集 的 网 格 化 弹 
P 性 缩放 到 ( 单 ) 多 机 系统 ， 充 分 发 挥 了 集群 的 规模 化 优势 。 
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射电 干涉 阵列 得 到 的 是 非 均匀 采样 的 可 见 度 数据 ,在 应 用 快速 传 里 时 变换 (Fast Fourier 
Transform, FFT) 对 可 见 度 数据 进行 成 像 前 ， 需 使 用 网 格 化 方法 将 实际 的 采样 数据 重 采样 到 
一 个 均匀 划分 的 网 格 上 。 当前 网 格 化 主要 是 使 用 基于 卷 积 的 网 格 化 方法 ,， 卷 积 网 格 化 过 程 的 
f 实质 是 矩阵 相 乘 ， 在 当 数 据 量 较 大 时 ， 网 格 化 计算 是 非常 耗 时 。 
:一 近年 来 , 天 文学 家 们 在 提升 可 见 度 数据 网 格 化 算法 性 能 做 了 很 多 研究 .其 中 W-projection 
< 算法 是 目前 广泛 使 用 的 网 格 化 方法 ， 由 于 该 算法 仅 校准 W 项 ， 并 没有 校准 方向 相关 效应 的 
A 项 , 当天 线 彼此 相距 较 远 , W 项 的 尺寸 可 能 会 变 得 很 大 , 使 该 算法 效率 低下 且 占 用 内 存 口 。 
通过 将 每 个 可 见 度 数据 的 w 值 投影 到 邻近 的 w 平面 的 W-Stacking 算法 ， 可 以 显著 提高 网 格 
化 性 能 ,但 是 需要 耗费 额外 的 内 存 站 。 如 果 考 虑 方向 相关 效应 ， 网 格 化 的 计算 难度 将 进一步 
增加 ， 同 时 修正 方向 相关 效应 4 项 和 W 项 被 称 为 AW-projection 网 格 化 算法 中 。 在 数值 分 析 
领域 ，Barnett 等 人 四 提出 基于 “半圆 指数 ” 卷 积 核 的 非 均 匀 侍 里 叶 变 换 库 (Non-uniform Fast 
Fourier Transform, NUFFT)， 将 FFT 推广 到 离散 化 的 网 格 数 据 中 。 首 次 将 NUFFT 应 用 到 射 
EE 天 文中 的 Nifty-gridder 算法 , 采用 共享 内 存 和 多 线程 技术 , 进一步 优化 W-Stacking 算法 。 
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综 上 所 述 , 网 格 化 算法 的 改进 和 细 化 都 需要 计算 更 多 的 卷 积 核 , 卷 积 计算 占据 网 格 化 算 
法 开销 的 主要 部 分 。 虽 然 采 用 多 核 CPU 和 GPUPI， 可 以 实现 并 行 计算 ， 提 高 算法 性 能 ， 但 
基于 Python 实现 的 上 述 网 格 化 方法 主要 局 限于 NumPy 多 维 数组 计算 , 难以 适应 数据 的 海量 
性 和 实时 处 理 需 求 。 近年 来 数组 Dask.Array 的 提出 , 为 超大 和 矩阵 的 数值 计算 开辟 了 新 途径 
Jamie Farnes 等 人 四 采用 Dask 并 行 框 架 外 ， 配 合 Pipeline 技术 ， 测 试 LOFAR 数据 集 ， 使 得 
原本 需要 11 个 小 时 才能 完成 整个 成 像 流 程 的 串 行 化 代码 ， 缩 短 至 8 分 钟 ， 大 大 减少 了 干涉 
成 像 所 需 的 时 间 。 本 文 提出 基于 Dask 并 行 加 速 的 射电 干涉 可 见 度数 据 卷 积 网 格 化 方法 ， 在 
并 行 计算 的 基础 上 兼顾 系统 的 弹性 缩放 ， 主 要 特点 是 以 Dask.Array 算 阵 分 块 存 储 和 计算 为 
核心 ， 封 装 Nifty-gridder 卷 积 网 格 化 算法 提供 的 Python 接口 ， 采 取 数 据 分 块 和 延迟 计算 ， 
提高 了 数值 计算 效率 ， 配 合 Dask 的 分 布 式 调度 策略 ， 实 现 了 网 格 化 算法 从 单机 到 集群 的 迁 
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其 中 (1，m， nn) 是 观测 方向 的 余弦 坐标 ，(wyo，wwo，wa ) 是 天 线 p 和 gq 组 成 的 基线 坐标 ， 


CS G 和 A 项 分 别 是 琼斯 矩阵 参数 化 的 方向 无 关 和 方向 相关 效应 。 在 小 场 近似 的 条 件 下 ， 指 数 
= 中 的 wa(n 一 1) 趋 近 于 零 ， 可 见 度 和 天 空 亮度 近似 为 二 维 傅 里 叶 变 换 关 系 。 由 于 基线 uv 轨迹 
的 不 规则 性 , 可 见 度数 据 并 非 等 间隔 离散 采样 ,直接 对 干涉 测量 方程 进行 侍 里 叶 反 演 的 计算 

代价 是 非常 昂贵 的 。 为 了 应 用 FFT 算法 成 像 ， 可 见 度数 据 必须 重新 采样 到 规则 化 的 笛 卡尔 

网 格 中 (Gridding)。 

在 成 像 流程 图 H, 不 同 的 光谱 频率 〈 即 图 像 通 道 ) 测量 所 得 的 可 见 度数 据 可 以 独立 处 
理 。 一 个 图 像 通道 通常 对 应 于 一 个 或 多 个 数据 通道 。 成 像 通常 从 空白 的 天 空 模型 开始 欠 代 ， 
经 过 网 格 化 和 傅 里 叶 逆 变 换 运 算 ， 一 个 或 多 个 明亮 的 源 可 能 会 掩盖 周围 微弱 的 光源 ， 使 用 
CLEAN 算法 提取 明亮 点 源 到 天 空 模型 中 。 与 网 格 化 相反 的 过 程 是 对 天 空 模型 进行 快速 傅 里 
叶 变换 ， 即 从 天 空 模型 计算 模型 的 可 见 度 ， 这 被 称 为 去 网 格 化 (Degridding)。 测 量 可 见 度 减 
去 模型 可 见 度数 据 是 为 了 进一步 提取 微弱 光源 .重复 网 格 化 和 去 网 格 化 , 直到 天 空 模型 收敛 。 
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图 1 射电 干涉 成 像 流程 


sky-image 


Fig.1 The imaging pipeline of radio interferometry 


干涉 成 像 是 射电 天 文中 的 关键 步骤 。 简 单 地 将 可 见 度 数据 插值 到 邻近 的 网 格 中 会 导致 严 
重 的 伪 影 ， 特 别 是 图 像 混 三 ( 视 场 外 的 强 光源 甚至 噪声 混 堆 到 视 场 中 )。 为 了 抑制 图 像 伪 影 的 
副作用 ， 通 常 采用 可 见 度数 据 与 网 格 化 函数 ( 卷 积 核 ) 进 行 卷 积 ， 然 后 再 重 采 样 到 网 格 中 ， 可 


以 提供 抗 重 倒 效果 。 由 于 卷 积 核 的 窗 函 数 特 性 ,使 得 边界 处 的 图 像 裁剪 误差 比 中 心 位 置 要 高 
出 几 个 数量 级 ， 产 生 了 较 大 的 脏 图 。 这 需要 脏 图 与 修正 函数 相 乘 来 抵消 卷 积 核 产 生 的 误差 ， 
从 而 获得 正确 的 光 通 量 ， 且 该 修正 函数 通常 是 卷 积 核 的 传 里 叶 逆 变换 的 倒数 。 相 比 于 


W-Stacking f 


1&, Nifty-gridder 为 提高 卷 积 核 的 计算 精度 做 了 以 下 改进 : (1) 沿 着 w 轴 ， 对 所 


有 的 可 见 度数 据 网 格 化 到 三 维 wvw 空 间 内 ,而 不 是 将 每 个 可 见 度 数据 的 w 值 投影 到 邻近 的 w 


平面 ; (2) 改 进 后 的 修正 函数 使 脏 图 的 FFT 和 DFT 之 间 的 差异 值 最 小 化 (DFT 是 无 损 变换 )， 
因此 获得 更 高 精度 的 脏 图 。 


3 卷 积 网 格 化 实现 


3.1 测量 集 的 并 行 读 取 和 分 块 


相 比 于 NumPyndarray 数组 ，Dask.Array 具有 如 下 优势 : (1) 支 持 将 超大 数组 分 割 成 许多 
个 NumPy.ndarray 小 数组 ; (2) 采 用 阻塞 算法 能 在 比 内 存 大 的 数组 上 支持 多 核 并 行 计 算 。 此 外 
我 们 利用 Xarray 实现 矩阵 的 一 致 性 分 块 (chunksize)， 相 关 字 段 的 数据 可 以 轻易 地 转化 为 


Dask.Array 类 


型 。 对 于 ( 单 ) 多 个 测量 集 文件 , 统一 将 路 径 信 息 放 入 列表 对 象 中 , 使 用 Dask.Bag 


分 布 式 加 载 , 然后 按照 测量 集中 的 FIELD_ID #11 DATA_DESC_ID 字段 分 组 ,实现 并 行 加 载 。 


在 本 实验 中 整个 数据 集 划 为 四 个 子 数据 集 : (0.0，0 1，0 2 和 0 _3)。 以 子 集 0_1 AB, F 


分 重要 字段 及 数据 类 型 如 下 表 1 所 示 。 


Dimensions: 


表 1 Xarray 数据 集 定义 的 部 分 相关 实验 数据 


Tab.1 Xarray dataset definitions for some related experiment data 


(ant: 27, chan: 64, corr: 2, row: 413696, uvw: 3, xyz: 3) 


Coordinates: 


ROWID 


Dimensions without coordinates: ant, chan, corr, row, UVW, XyZ 


Data variables: 


ANTENNAIL 


ANTENNA2 


FLAG 


DATA 


UVW 


WEIGHT 


CHAN_FREQ 


Attributes: 


FIELD_ID: 


DATA_DESC_ID: 


(row) int32 dask.array<chunksize=(20000,), meta=np.ndarray> 


(row) int32 dask.array<chunksize=(20000,), meta=np.ndarray> 


(row) int32 dask.array<chunksize=(20000,), meta=np.ndarray> 


(row, chan, corr) bool dask.array<chunksize=(20000, 64,2),meta=np.ndarray> 


(row, chan, corr) complex64 dask.array<chunksize=(20000, 64, 2), meta=np.ndarray> 


0 


1 


3.2 ”网 格 化 方法 的 并 行 实现 


(row, uvw) float64 dask.array<chunksize=(20000, 3), meta=np.ndarray> 


(row, corr) float32 dask.array<chunksize=(20000, 2), meta=np.ndarray> 


(chan) float64 dask.array<chunksize=(64,), meta=np.ndarray> 


分 布 式 计算 是 解决 海量 数据 的 有 效 途径 ，Dask 并 行 计算 框架 提供 了 灵活 多 变 的 分 布 式 


调度 方式 。 由 于 Dask 任务 调度 方式 和 | 


便 可 以 使 算法 在 


“—~ 


单 ) 多 机 以 多 ( 线 ) 进 程 的 方式 弹性 


多 个 MS 文人 


源 ( 即 多 个 Sub-dataset)， 基 于 子 数 据 集 的 任务 调度 更 进 


F 的 物理 性 分 离 有 利于 使 用 多 进 


化 的 并 行 流程 ,使 用 高 阶 函 数 Dask. Array.blockwise 圭 


护 


TP Ee ASIA ) 


LEIE o FES A Wel BE GOEL 
程 并 行 读 取 数据 集 。MS 文件 通常 
步 细 粒度 化 整个 Nifty-gridder 网 格 


| 装 和 调用 Nifty-gridder 的 Python 接口 ， 


实现 了 基于 子 块 的 并 行 计算 以 及 协调 子 块 的 缩聚 和 拼接 操作 (图 
进程 之 间 的 传输 成 本 ， 数 值 计 算 采 用 多 线程 的 执行 方式 计算 脏 图 。Nifty-gridder 算法 的 执行 


过 程 如 下 : 


j 户 只 需 切换 调度 方式 ， 
展 ， 但 需要 根据 算法 的 特点 ， 选 择 合理 
的 任务 调度 方式 ， 以 获取 最 佳 的 计算 性 能 。 本 文 使 用 最 为 复杂 的 dask.distributed 调度 方式 在 
两 台 机 器 节点 执行 Nifty-gridder 网 格 


2.a) 采 用 多 进 


2.b)。 为 避免 


程 的 执行 方式 : 
包含 多 个 射电 


Dask.Array 在 


C1) 沿 着 w 轴 确 定 Nw 个 采样 平面 ， 并 均匀 分 布 到 w 轴 (从 wo~wy,-_1); 


(2) 沿 着 w 轴 


将 可 见 度数 据 网 格 化 到 w 平 


(3) 初始 化 Ne * My 的 零 矩 阵 1， 对 每 一 个 w = wi 平面 有 : 


a) 
b) 


c) 


将 每 个 w 平面 再 进行 uv 网 格 化 ， 然 后 执行 二 维 


传 


叶 逆 变换 ; 


裁剪 掉 FFT 图 像 的 外 半 部 ， 然 后 乘 上 e2miwiC1 一 


将 上 述 结果 累加 到 矩阵 I 中 ， 
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(4) 修正 函数 乘 以 矩阵 1， 得 到 最 终 的 脏 图 。 


Sub-dataset 
chunk chunk 
Receive MS Process Receive MS dataset dataset 


on dask-worker = Space on dask-worker | Thread Space 
gridding gridding 
kernel 3 kernel 
Sub-dataset ...... LoadMSand 
| Split by Source Sub-dataset 
| | iFFT 


Gridding 


S TESO ea Gridding Process on Gridding 
z IFFT each Sub-dataset ~~ and 
iFFT 


Array of ditry 
image 


ask Array ...... List of Dirty Image Array ...... Dask.Array 


(a) Task scheduling on dask.distributed (b) Parallel computation in chunk 


图 2 分 布 式 任务 调度 和 Nifty-gridder 网 格 化 算法 


Fig. 2 Distributed task scheduling and the Nifty-gridder algorithm 


4 实验 结果 和 分 析 


4.1 实验 环境 


实验 的 数据 集 ! 来 源 于 2010 年 8 月 23 A, 由 甚大 型 Karl G. Jansky 干涉 阵列 对 超新星 遗 
WE G055.7+3.4 进行 长 达 8 小 时 的 观测 。 该 阵列 采用 D- 型 配置 ， 观 测 频 率 范 围 为 1G-2GHz， 
履 盖 所 有 可 用 的 工 -波段 。 实 验 的 硬件 环境 为 两 台 高 性 能 服务 器 : Intel Xeon CPU E5-2660 v4 
CPU @ 3.4GHz 处 理 器 (56 %), 512GB RAM。 使 用 Common Astronomy Software 
Applications(CASA 5.6.2) 进 行 数据 结果 的 验证 。 


4.2 Dask 并 行 加 速 和 实验 结果 


以 四 个 子 数 据 集 为 例 ，chunksize 设置 为 20000 行 ， 经 网 格 化 处 理 生成 脏 图， 使 用 可 见 
度数 据 的 行 数 度量 数据 集 的 体积 ， 在 同一 软 硬 件 环境 下 比较 Dask.Array 和 NumPy 版 本 
Nifty-gridder 算法 的 运行 时 间 ( 单 位 ， 秒 )， 实 验 结果 如 下 : 


表 2 Nifty-gridder 网 格 化 执行 时 间 的 比较 (Dask.Array vs. NumPy) 


http://casa.nrao.edu/Data/EVLA/SNRGS55/SNR_G55_10s.calib.tar.gz 
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Tab.2 The execution times of the Nifty-gridder compared Dask.Array to NumPy 


Execution Time (Dask.Array) 
Sub-dataset Volume(row) Execution Time (NumPy) Speedup Ratio 
CPU Time(s) Wall Time(s) 
0_0 39274 2.8 2.68 1.57 0.59 
0_1 413696 38.5 4.95 16.18 3.27 
0_2 412696 35.9 4.86 14.26 2.93 
0_3 414974 39.4 4.95 16.16 3.26 


YE: Speedup Ratio= Execution Time (NumPy) / Wall Time(s), H. Wall Time 为 程序 的 实际 执行 时 间 。 

从 表 2 可 知 ,基于 Dask. Array 改进 的 Nifty-gridder 算法 ,其 CPU Time LYK F Wall Time, 
说 明 对 于 计算 密集 型 问题 ， 使 用 多 核 计 算 并 行 效果 显著 ， 明 显 降 低 程 序 的 运行 时 间 。 以 0_0 
和 0_1 数据 集 的 对 比分 析 为 例 : 即使 将 可 见 度 数据 体积 增 大 10.5 倍 (<*413696/39274)， 相 应 
的 执行 时 间 Wall Time 仅仅 增加 1.85 f(~4.95/2.68), 且 加 速 比 进一步 提高 。 然 而 Dask.Array 
是 在 NumPy 的 基础 上 增加 了 一 层 复杂 的 设计 ， 对 于 较 小 的 数据 体积 (0_1 数据 集约 占用 
40MB), NumPy 可 能 是 正确 的 选择 ， 相 反 ， 这 恰恰 说 明了 Dask.Array 适宜 处 理 超大 型 矩阵 
的 数值 计算 。 

Dask 允许 跨 集群 提交 Python 函数 以 实现 基于 任务 的 并 行 ， 从 而 生成 大 量 可 以 监视 、 控 
制 和 计算 的 Future 对 象 。 对 于 复杂 的 程序 处 理 流程 ， 动 态 的 可 视 化 监控 有 助 于 了 解 算 法 的 
性 能 瓶颈 ， 实 验 执行 过 程 中 的 实时 性 能 监控 如 图 3 所 示 。 
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图 3 网 格 化 流程 中 任务 流 的 实时 状态 


Fig. 3 The real-time state of task stream in the gridding process 
为 了 说 明 Dask 并 行 框架 的 优越 性 ， 通 过 增加 测量 集 的 输入 量 和 限定 每 台 机 器 内 存 占 有 
量 并 确保 实验 环境 一 致 。 从 系统 资源 利用 率 角 度 分 析 并 比较 基于 Dask.Array 和 NumPy 的 
Nifty-gridder 算法 性 能 。 由 图 4 可 知 ， 无 论 是 资源 利用 率 的 峰值 和 平均 值 ， 相 比 于 NumPy 
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版 本 ，Dask.Array 类 型 的 网 格 化 算法 明显 占有 更 低 CPU 利用 率 和 内 存 占 有 率 ， 但 却 能 获得 
更 快 的 网 格 化 执行 时 间 ( 见 表 2)。 主 要 是 因为 Dask.Array 数组 是 采取 分 块 加 载 和 延迟 计算 ， 
尚 不 具备 计算 条 件 的 子 块 会 驻 留 磁盘 ， 以 节约 系统 资源 ， 而 满足 计算 条 件 的 子 块 则 被 送 入 内 


存 并 行 执行 ， 相 反 NumPy 数组 必须 全 部 加 载 到 内 存 ， 导 致 较 高 内 存 的 持 有 率 。 
CPU 
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图 4 网 格 化 流程 中 CPU 和 内 存 的 使 用 情况 对 比 (Dask.Array vs. NumPy) 


Fig.4 CPU and memory usage in the gridding process compared Dask.Array to NumPy 

为 了 进一步 验证 代码 的 正确 性 ， 使 用 标准 的 CASA 软件 对 该 数据 集 进行 成 像 ， 生 成 的 
脏 图 (图 5 左 ) 与 实验 结果 (图 5 右 ) 进 行 对 比 ， 两 幅 灰 度 图 中 的 灰白 色 点 代表 观测 源 ， 可 以 发 
现 正确 识别 出 射电 源 的 分 布 位 置 。 


> https://casaguides.nrao.edu/index.php?title=VLA_CASA_Imaging-CASAS.7.0 
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图 5 CASA 和 实验 结果 的 脏 图 对 比 
Fig.5 Comparison of dirty image from CASA and experimental result 
5 总 结 


高 性 能 分 布 式 并 行 计算 已 成 为 射电 干涉 成 像 过 程 中 应 对 高 分 汰 率 和 大 视 场 干涉 阵列 产 
生 的 海量 数据 的 必要 方法 。 可 见 度数 据 的 网 格 化 和 去 网 格 化 是 成 像 的 重要 组 成 部 分 , 网 格 化 
并 行 加 速 无 颖 对 于 提高 整个 成 像 过 程 速度 有 重要 意义 。 本 文 使 用 了 开源 的 Dask 分 布 式 计算 
框架 结合 Nifty-gridder 实现 了 测量 集 的 分 布 式 加 载 和 并 行 网 格 化 加 速 过 程 , 充分 发 挥 了 集群 
的 规模 化 优势 ， 提 高 了 多 核 CPU 利用 率 。 干 涉 成 像 过 程 中 包含 多 个 复杂 的 处 理 流程 ， 都 涉 
及 矩阵 的 数值 计算 ， 而 Dask.Array 可 以 高 效 地 处 理 多 维 超大 和 矩阵 的 数值 计算 ， 因 此 下 一 步 
的 工作 考虑 基于 Dask 实现 去 网 格 化 、 校 准 、 成 像 等 算法 的 并 行 加 速 。 
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A distributed gridding implementation method for radio 


interferometric visibilities based on DASK 


Li Shanshan', Luo Kaida', Wei Shoulin'’, Dai Wei'”, Liang Bo!” 
(1, Faculty of Information Engineering and Automation, Kunming University of Science and Technology, Kunming 650500, China; 
2, Key Laboratory of Applications of Computer Technology of the Yunnan Province, Kunming 650500, China) 

Abstract: Fast Fourier Transform (FFT) has better performance than Discrete Fourier 
Transform, which is the fundamental imaging algorithm of radio interferometry. However, 
because of the irregular sampling of antenna array, it is necessary to use gridding algorithms to 
resample visibilities to regular grids, so that FFT can be applied. The convolutional gridding in 
radio interferometric imaging is characterized by intensive and iterative computations. Especially 
in the case of massive visibility data processing, high-performance gridding computing is 
particularly important to accelerate the whole imaging process. In order to alleviate the pressure of 
data processing, the DASK parallel computing framework is extended and applied on the existing 
gridding algorithm which supports multi-core parallelism but processes whole blocks of data. Not 
only can the data be partitioned and distributed to multiple threads to improve the efficiency of 
numerical computation, but also the dynamic distributed task scheduling strategy can optimize the 
real-time workflow of gridding. The experimental results show that the multi-core utilization rate 
is significantly improved and the performance of gridding algorithm can be further enhanced even 
if the volume of visibility is increased. Distributed task scheduling can flexibly scale the gridding 
task of (single) multi-measurement set to (single) multi-machine system, which gives full play to 
the scale advantage of clustering. 
Key words: Gridding; Interferometric imaging; Distributed Computing; Parallel Computing; 


Dask 


